{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from [Derivative Approximation by Finite Differences](https://www.geometrictools.com/Documentation/FiniteDifferences.pdf) by David Eberly,  additional text and SymPy examples by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalization of Taylor FD operators\n",
    "\n",
    "In the last lesson, we learned how to derive a high order FD approximation for the second derivative using Taylor series expansion. In the next step we derive a general equation to compute FD operators, where I use a detailed derivation based on [\"Derivative Approximation by Finite Differences\" by David Eberly](https://www.geometrictools.com/Documentation/FiniteDifferences.pdf)\n",
    "\n",
    "## Estimation of arbitrary FD operators by Taylor series expansion\n",
    "\n",
    "We can approximate the $d-th$ order derivative of a function $f(x)$ with an order of error $p>0$ by a general finite-difference approximation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) = \\sum_{i=i_{min}}^{i_{max}} C_i f(x+ih) + \\cal{O}(h^{d+p})\n",
    "\\end{equation}\n",
    "\n",
    "where h is an equidistant grid point distance. By choosing the extreme indices $i_{min}$ and $i_{max}$, you can define forward, backward or central operators. The accuracy of the FD operator is defined by it's length and therefore also the number of \n",
    "weighting coefficients $C_i$ incorporated in the approximation. $\\cal{O}(h^{d+p})$ terms are negelected. \n",
    "\n",
    "Formally, we can approximate $f(x+ih)$ by a Taylor series expansion:\n",
    "\n",
    "\\begin{equation}\n",
    "f(x+ih) = \\sum_{n=0}^{\\infty} i^n \\frac{h^n}{n!}f^{(n)}(x)\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Inserting into eq.(1) yields\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) &= \\sum_{i=i_{min}}^{i_{max}} C_i \\sum_{n=0}^{\\infty} i^n \\frac{h^n}{n!}f^{(n)}(x) + \\cal{O}(h^{d+p})\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "We can move the second sum on the RHS to the front\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) &= \\sum_{n=0}^{\\infty} \\left(\\sum_{i=i_{min}}^{i_{max}} i^n C_i\\right) \\frac{h^n}{n!}f^{(n)}(x) + \\cal{O}(h^{d+p})\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "In the FD approximation we only expand the Taylor series up to the term $n=(d+p)-1$\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) &= \\sum_{n=0}^{(d+p)-1} \\left(\\sum_{i=i_{min}}^{i_{max}} i^n C_i\\right) \\frac{h^n}{n!}f^{(n)}(x) + \\cal{O}(h^{d+p})\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "and neglect the $\\cal{O}(h^{d+p})$ terms\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) &= \\sum_{n=0}^{(d+p)-1} \\left(\\sum_{i=i_{min}}^{i_{max}} i^n C_i\\right) \\frac{h^n}{n!}f^{(n)}(x)\\\\\n",
    "\\end{align}\n",
    "\n",
    "Multiplying by $\\frac{d!}{h^d}$ leads to the desired approximation for the $d-th$ derivative of the function f(x):\n",
    "\n",
    "\\begin{align}\n",
    "f^{(d)}(x) &= \\frac{d!}{h^d}\\sum_{n=0}^{(d+p)-1} \\left(\\sum_{i=i_{min}}^{i_{max}} i^n C_i\\right) \\frac{h^n}{n!}f^{(n)}(x)\\\\\n",
    "\\end{align}\n",
    "\n",
    "Treating the approximation in eq.(2) as an equality, the only term in the sum on the right-hand side of the approximation that contains $\\frac{h^d}{d!}f^{d}(x)$ occurs when $n = d$, so the coefficient of that term must be 1. The other terms must vanish for there to be equality, so the coefficients of those terms must be 0; therefore, it is necessary that\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i=i_{min}}^{i_{max}} i^n C_i=\n",
    "\\begin{cases}\n",
    "0, ~~ 0 \\le n \\le (d+p)-1 ~ \\text{and} ~ n \\ne d\\\\\n",
    "1, ~~ n = d\n",
    "\\end{cases}\\nonumber\\\\\n",
    "\\end{equation}\n",
    "\n",
    "This is a set of $d + p$ linear equations in $i_{max} − i_{min} + 1$ unknowns. If we constrain the number of unknowns to be $d+p$, the linear system has a unique solution. \n",
    "\n",
    "- A **forward difference approximation** occurs if we set $i_{min} = 0$\n",
    "and $i_{max} = d + p − 1$. \n",
    "\n",
    "- A **backward difference approximation** can be implemented by setting $i_{max} = 0$ and $i_{min} = −(d + p − 1)$.\n",
    "\n",
    "- A **centered difference approximation** occurs if we set $i_{max} = −i_{min} = (d + p − 1)/2$ where it appears that $d + p$ is necessarily an odd number. As it turns out, $p$ can be chosen to be even regardless of the parity of $d$ and $i_{max} = (d + p − 1)/2$.\n",
    "\n",
    "We could either implement the resulting linear system as matrix equation as in the previous lesson, or simply use a `SymPy` function which gives us the FD operators right away."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import SymPy libraries\n",
    "from sympy import symbols, differentiate_finite, Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1st order forward operator 1st derivative:\n",
      "(-f(x) + f(h + x))/h\n",
      " \n",
      "1st order backward operator 1st derivative:\n",
      "(f(x) - f(-h + x))/h\n",
      " \n",
      "2nd order center operator 1st derivative:\n",
      "(-f(-h + x) + f(h + x))/(2*h)\n",
      " \n",
      "2nd order center operator 2nd derivative:\n",
      "(-2*f(x) + f(-h + x) + f(h + x))/h**2\n",
      " \n",
      "4th order center operator 2nd derivative:\n",
      "(-30*f(x) - f(-2*h + x) + 16*f(-h + x) + 16*f(h + x) - f(2*h + x))/(12*h**2)\n",
      " \n"
     ]
    }
   ],
   "source": [
    "# Define symbols\n",
    "x, h = symbols('x h')\n",
    "f = Function('f')\n",
    "\n",
    "# 1st order forward operator for 1st derivative\n",
    "forward_1st_fx = differentiate_finite(f(x), x, points=[x+h, x]).simplify()\n",
    "print(\"1st order forward operator 1st derivative:\")\n",
    "print(forward_1st_fx)\n",
    "print(\" \")\n",
    "\n",
    "# 1st order backward operator for 1st derivative\n",
    "backward_1st_fx = differentiate_finite(f(x), x, points=[x, x-h]).simplify()\n",
    "print(\"1st order backward operator 1st derivative:\")\n",
    "print(backward_1st_fx)\n",
    "print(\" \")\n",
    "\n",
    "# 2nd order centered operator for 1st derivative\n",
    "center_1st_fx = differentiate_finite(f(x), x, points=[x+h, x-h]).simplify()\n",
    "print(\"2nd order center operator 1st derivative:\")\n",
    "print(center_1st_fx)\n",
    "print(\" \")\n",
    "\n",
    "# 2nd order FD operator for 2nd derivative\n",
    "center_2nd_fxx = differentiate_finite(f(x), x, 2, points=[x+h, x, x-h]).simplify()\n",
    "print(\"2nd order center operator 2nd derivative:\")\n",
    "print(center_2nd_fxx)\n",
    "print(\" \")\n",
    "\n",
    "# 4th order FD operator for 2nd derivative\n",
    "center_4th_fxx = differentiate_finite(f(x), x, 2, points=[x+2*h, x+h, x, x-h, x-2*h]).simplify()\n",
    "print(\"4th order center operator 2nd derivative:\")\n",
    "print(center_4th_fxx)\n",
    "print(\" \")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, the underlying algorithm also supports variable grid spacings, because it is not based on Taylor series expansion but Lagrange polynomials. For more details, I refer to the paper [\"Calculation of weights in finite difference formulas\" by Bengt Fornberg](https://amath.colorado.edu/faculty/fornberg/Docs/sirev_cl.pdf)\n",
    "\n",
    "An alternative to using `SymPy` for the calculation of FD operator weights is the **DEVITO** package:\n",
    "\n",
    "[\"https://www.devitoproject.org/\"](https://www.devitoproject.org/)\n",
    "\n",
    "It does not only calculate FD stencils, but also automatically optimizes the performance for the given hardware."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* How to compute Finite-Difference operators of arbritary derivative and error order\n",
    "* Symbolic computation of FD operators with `SymPy`"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
